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TECHNICAL MEMORANDUM X-53997 


SOME APPLICATIONS AND LIMITATIONS 
OF THE FAST FOURIER TRANSFORM 

INTRODUCTION 


The analysis of a time series is very important in many experiments 
(such as vibrational testing of structures). Seemingly, the two functions 
desired most often are the power spectral density function and the correlation 
function. Theoretically, if one is known exactly, so is the other. In 
practice, small errors in one may cause large errors in the other. This 
discussion includes power spectral density functions, autocorrelation 
functions, cross-spectral density functions , cross-correlation functions, 
and other related functions. 

There are three separate aspects of the spectral analysis, 

(1) definition of spectrum, (2) calculation of the defined spectrum, and 
(3) interpretation of the calculated spectrum. All three problems are 
considered here. 

The most direct route from a time series (discrete or continuous) 
to a spectral density function is the Fourier transform. In the past, the 
computer time required to calculate a Fourier transform made this approach 
impractical; as a result, most approaches were indirect. The correlation 
function for a limited number of lags was calculated, then the Fourier 
transform of the correlation function was calculated. Recent developments 
in procedures for complex arithmetic have made the calculation of Fourier 
transforms of discrete time series feasible. Because older methods have 
been in use longer, change to the direct Fourier transform method is slow. 

This discussion explains the uses and limitations of the new approach 
using the Fast Fourier Transform. 



FOURIER TRANSFORM 


A detailed discussion of the Fourier transform, with the necessary 
and sufficient conditions for the existence of the transforms, is found in 
many textbooks. Here, the existence of the integrals or sums will be 
assumed to exist. In the Appendix, the procedure known as the Fast 
Fourier Transform is explained more fully. 

If x(t) is a time series and if the two integrals 

CO 

X(f) = / x(t)e - ^ 27r ft dt U) 

-OO 

and 

°o 

x (t) = f X(f)e i2n ft df (2) 


exist, X(f) is called the Fourier transform of x(t) , and x(t) is called 
the inverse Fourier transform of X(f) . They, X(f) and x(t) , are 
frequently called transform pairs. 

In general, the time series that are recorded for analysis are finite 
and the existence of the preceding integrals are guaranteed by the assump- 
tion that the series is zero outside some time interval (usually 0 ^ t < T). 
Under these assumptions these integrals become 

X(f) = / x(t)e' j27rft dt (3) 

0 


x(t) = / X(f)e j27r ft df , (4) 

-F 

where F is an upper bound of the frequencies found in x(t) . Also, in 
practice, the time series is a discrete time series (or must be made into 
one) of N values usually equally spaced. Therefore, the above integrals 
become 


2 



(5) 


N-i -327 r ft. 

X(f) = (At) Yj x(t.) e 1 

i=0 1 


with 


T = N • At 


and 


x(t ) = (Af) 
k 


N-i 


Z 


i=0 


X(f.) 


+j2,r f.t k 

e 


f. = i(Af) 


(6) 


where the last expression is summed on frequency. These sums are 
probably more familiar as 

X(f) = (At) Y x.(t.) ^cos 27r ft^ - j sin 27r ft. J (7) 

and 

x(t) = (Af) Y X(l) ^cos 27T Tt + j sin 27r f.t J (8) 

In the first formula, one may calculate X(f) for any f by simply inserting 
that value of f in the sum and summing the values. In the second equation, 
the same is true of any t value if X(f) is known for all f . These 
transforms are usually referred to. as Discrete Fourier Transforms, or 
(DFT). If the sampling rate is sufficient, the value of the Discrete Fourier 
Transform should agree very closely with the Fourier transform evaluated 
at the same f value. 

The Fast Fourier Transform (FFT) is just the DFT with some added 
restrictions and some innovations in the method of calculation. These 
restrictions are (1) N (the number of points) must be a composite number; 

(2) the calculated values of f are 0, ^ , . . . , " ; and (3) the 
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points must be equally spaced. For the limited values where X(f) is 
calculated, the FFT and the DFT should agree to the accuracy of the program 

cl 

used. A great many of the programs available insist that N = 2 for some 
integer a . It seems to be desirable in all cases that N have no prime 
factor greater than 5. 

Since the transform is the integral part of this study, it must be 
calculated with a reasonable degree of accuracy. However, the function of 
concern is the power spectrum and the discussion of accuracy will be 
included in that section. 


POWER SPECTRUM 

There are several definitions of power spectrum in the literature [ 1 ] . 
In theory, the direct definition of the power spectral density function of x(t) 
is 


S x (f) = lim ± [x T (f)X T (f)] , (9) 

where 

X (f) = / x(t) e _;i27rft dt (10) 

-T 

and X T (f) is the complex conjugate of X T (f) . Since x(t) is a finite 
series in the cases considered here, the definitions will be 

S X (f) = T [^ (f >X(f)] (Ha) 

and 

T 

X(f) = f X(t) e"^ 27r ft dt (lib) 

0 
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This definition is for a time series where t runs from 0 to T and f runs 
N N 

from - — Af to — Af in the usual Discrete Fourier Transform approach, 

Z Z 

but f runs from 0 to N (Af ) in the Fast Fourier Transform approach. 

Both of these are effectively two-sided spectral functions. It can be shown 

that if X(f) is calculated from - ^ Af to ~ Af , X(kAf) =X(-kAf) 

Z Z 

and if f is calculated at 0, Af NAf , then X(kAf) =X[(N-k)Af]. 

Therefore in the first case, S (f) is symmetric about f = 0 , and in the 

N 

second case it is symmetric about — (Af) . 

z 

Some writers call the function of equations (lia) and (lib) the 
periodogram of x(t) [2] and take m time records of length T and average 
these to get the spectral density function. Their definition would be 


S (f) 
x 


m 


m 

2 V'> 

1=1 l 


(12) 


where S (f) is the spectral density function of the ith time series by the 

X . 

1 

previous definition or the periodogram by their definition. In actual 
problems, it is common practice to have a one-sided spectral function, 
defined as 


G (f) = 2S (f) , 0 < f < ^ (Af) 

xx 2 


(13) 


In this report, the spectral density function will be as defined by 
equations (11a) and (lib). If some modification is desired, it will be 
called the modified spectral function. 

There are many ways to modify the spectral density function. One is 
given above where m series of length T are used and their spectral 
density functions are averaged. The more common modifications are to 
weight the time series before finding the Fourier transform, or by weighting 
the calculated Fourier transform before calculating the spectral density 
function. Usually, these procedures are referred to as weighting, 
smoothing or filtering. Some of the common weighting functions or filters 
are given later in this report [1, 2, 3]. However, for deterministic 
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functions, the results are more accurate without the filtering when the time 
T and the sampling rate are correctly chosen. 

Some other useful functions, easily calculated from the FFT and 
power spectral density functions, are discussed here. If x(t) and y(t) 
are two time series of the same length and sample rate, the cross power 
spectral density function is 


S (f) = lim ~ {x(f) Y (f ) } 

^ X 00 

or for our approximation 

S (f) =^{X(f)Y(f)> 
xy T 


(14) 


(15) 


Also, the autocorrelation function and cross-correlation function are 
Fourier transforms of the power spectral density function and cross spectral 
density function, respectively, i.e. , 


R ( t ) = f S (f) e^ 27r fT df 
x J x 


and 


R 

xy 


(t) 


f S (f) e 327r fT df 

J xv 


(16) 


(17) 


Under the restrictions of the Appendix, if S^(f) is entered instead of 

the time series in the FFT routine, the autocorrelation function is the output; 
if S x y(f) is used for input, the cross-correlation function is the output. 

Another function frequently desired is the phase angle between two 
time series. If the cross power spectral density function S (f) is written 

as a complex function 


S (f) = C (f) - JQ (f) 
xy xy xy 


(18) 
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then the phase angle 6 is 


6 


tan 


. Q (f) 

-1 xy 

C (f) 
xy 


( 19 ) 


If the correlation functions are desired, it is recommended by some [i, 2] 
that the time series x(t) or y(t) be altered by addition of the same 
number of zeros as the number of values of the correlation function desired. 
For example, if the time series consists of 1024 points and 100 values of 
the correlation function are desired, the last 100 values would be replaced 
by zeros. With the altered time series, the correlation functions computed 
via the FFT yields better agreement to the direct correlation method. 


FILTERS OR WEIGHTS 


Two commonly used programs for vibration analysis are RAVAN [4] 
and STAN 1 . These use the Blackman- Tuckey method. 

In the Blackman-Tuckey method of calculating the power spectra, 
filtering is necessary even for periodic deterministic data. With the FFT 
approach, filters are not required for periodic data when the length of the 
record is chosen correctly. 

The filter or weight should be designed for the particular problem. 
However, some of the common filters are discussed here. 

Defining the time series x(t) to be zero outside a finite interval 
- T ^ t < T is equivalent to using a weight function w(t) =1 on this interval 
and zero elsewhere. Now, according to Fourier transform theory, the 
Fourier transform of a product w(t) x(t) is the convolution of W(f) with 
X(f). Mathematically, this can be expressed as the transform of w(f) x(t) 


oo 

/ X(f 0 ) W(f-f 0 ) df 0 

— OO 


( 20 ) 


The effect of using a short sample of the record is shown in Figure 1 [5] . 


i. Newberry, M. H.: The Statistical Analyzer Program (STAN). Internal 
Note IN-COMP-67-1, George C. Marshall Space Flight Center, September 
1967. 
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FUNCTIONS 


TRANSFORMS 



NOTE : Ignoring the early and late parts of the function has 

the effect of "blurring" its spectrum. As shown 
here, the broad peaks AA of the transform are little 
affected; the rapid oscillations BB almost disappear; 
and the narrow peak at zero frequency is reduced. 

This conclusion is physically reasonable since cutting, 
off the ends of the function (a) leaves most of the high 
frequency oscillations and simultaneously deletes the 
slow ones and halves the area under the curve. The 
function has been chosen symmetrically to avoid the 
need to represent complex values in the transform. 

Figure 1. Example of the convolution theorem. 

If the first line (a) is x(t) and its transform X(f) , and the second 
line (b) is the weight function w(t) and its transform W(f) , then the last 
line (c) is w(t) x(t) and its transform. Notice that the higher frequency 
components of the transforms of x(t) and w(t) x(t) are very similar. 
However, in this case, the lower frequency components are greatly altered. 
Also, the value at f = 0 is much lower in the weighted case. The process of 
convolution causes a cancellation when there is a rapid oscillation, as in 
this case. 
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For the Discrete Fourier Transform or Fast Fourier Transform, 
this same situation applies unless the frequencies are discrete and the time 
interval is chosen so that the calculated values are the actual frequencies 
present. In this case with deterministic functions, the results are very 
accurate. 

This problem of using a finite record cannot be avoided. Comparison 
of this analysis with that of a longer segment of the time series will often 
help one to see the actual situation since the narrower bandwidth tends to 
separate the peaks. 

One of the weights recommended is [ 2 ] 


w(t) 

w(t) 

w(t) 


I ( 1 -° os rrr) 


for 0 ^ t < 0. IT 


1 for 0. IT < t ^ 0. 9T 


1 

2 


1 + COS 7T 


(t - 0. 9T) 
0. IT 


for 0. 9T < t < T 


(21) 


If one considers the above weight function and compares its spectra or 
transform with the square weighting function above, one sees that the side 
lobes are much smaller. This means that the frequencies in these side 
bands contribute less to the value calculated for a given frequency [1]. 

Another weight function used [6] is 


w(tj) = sin 2 -j-jj (22) 

or x. is weighted by w(t.) . One of the advantages of this filter is that the 
Fourier coefficients a^'s and b^'s can be calculated without the weight 
function, and then the coefficients can be weighted so that 


a ' = ~ 
k 2 


~ a. , 2a, - a, . 

k-1 k k+1 


) 


(23) 
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and 


b k -i (- b k -i + 2 \-\ + i) < 24 

The document cited in Reference 6 recommends further smoothing by 
calculating 

p k =1 ( a k 2+b k 2 ) ■ (25 

where P is the average power due to the kth frequency component, and 

K. 

then calculating 


p k =i( p k -i + 2p k + p k + i ) (26) 

when noise is present. This last step is equivalent to weighting the auto- 

o ITT 

covariance function with cos 2 — in the STAN program. 

There are many other documented weighting or filtering approaches 
[ 2 , 5 , 6 ]. 

It should be realized that smoothing, filtering or weighting basically 
degrades the resolution. Sharp peaks that might occur in a particular record 
are smoothed out so that the analysis of different records of the same 
general happening are more nearly the same. The improvement is in 
stability or reproduciblity. The user should decide first whether he wants 
resolution or stability. Then, if he desires stability, the proper filter 
should be chosen for his problem. 

One of the most commonly used class of weights is a generalized 
class of cosine tapers. The transform of the generalized form is 
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W(2irf, m) 



where m = fT . This weight is capable [3] of reducing the contribution 

from adjacent frequencies separated by more than f = tjt • However, 

optimum choice of m requires the advance knowledge of the minimum 
separation of any two frequencies of interest. 


ACCURACY OF FFT 


For the FFT problem, the main sources of error are sampling rate, 
numerical integration, quantization and the number of calculations required. 

The numerical integration scheme approximates the Fourier integral 
by a series 

T 

f x(t) (cos 27T ft - j sin 2ir ft) dt = (At) 

0 

(28) 

where N • At = T and 

y. = x(t.) I" cos 2 tt ft. - j sin 2ir ft.l (29) 

l l L l lj 


’ y 0 + V" +y N-/ 


This particular form does not lend itself to easy theoretical analysis. How- 
ever, it is very similar to the trapezoidal integration formula 


(At) 



+ y. 


* + 


y N-l 



(30) 
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for which the error estimates are known. The upper bound to the error in 
this formula is [7] 


E 


max 


Ny" (At) 3 
12 


(31) 


where y" is the second derivative of x(t) [cos 2ir ft - j sin 2ir ft]. If it is 
assumed that x(t) can be expanded into a trigonometric series, the four 
terms involved for each 1, are (cos 2w f t) (cos 2ir ft) , 

- j(sin 27r f.t) (sin 27rft), -j(cos 27rf.t) sin 2w ft, sin 2ir f .t cos 2w ft . 

These terms are very similar and only one is treated here. Since the 
integrals of these terms vanish unless f = f. , consider y(t) = cos 2 2w ft. 

Then, y'(t) is -4tt f cos 2ir ft(sin 2w ft) , y"(t) = 8ir 2 f 2 (sin 2 27r ft - 
cos 2 27T ft) , and |y"(t)| ^ 871 2 ! 2 . Now, let T = 1 unit of time. Then, 

if k samples are taken on each cycle of y(t), N = kf and At = rjr . 

KI 

Hence, 

(kt) <8rt»)(^)' _ 

0.12 12k 2 


The integral of cos 2 2tt ft between 0 and 1 is ^ . Hence, the maximum 
percentage error is or roughly jj- . This indicates the maximum 


k 2 


error for 10 samples per cycle of the given frequency is 14 percent. In 
deterministic cases , the error is much smaller than the theoretical limit 
when the sample rate exceeds 5 samples per cycle. For lower rates, it is 
very close. In one problem with 3 samples per cycle, the error was 100 
percent. Thus, errors resulting from sample rate and errors resulting 
from numerical integration procedure are closely associated. 


The quantization errors must be treated statistically. If it is 
assumed that the quantization errors satisfy certain common statistical 
models, the errors in the final answers are very small. Two examples 
that will be discussed later indicate that this conclusion is valid in these 
cases. 
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The propagation and round-off errors have been treated by Welch [8] 
and have not been treated here. Using his assumptions, the percentage 
error was very low, usually less than 0. 1 percent. 

While error estimates are very worthwhile, experiments with 
functions where the answers are known or can be compared with some 
accepted standard are very desirable. 

The first case chosen was the deterministic function 
X(f) = 10 cos 20. 5 7rt , (f = 10.25) . The objective in this experiment was 
two-fold: (1) to show the effect of the record length T on resolution, 

(2) to show the accuracy of the method. 

The results when the time T was ~ , t > ^ > 1 , 2, and 4 seconds are 

8 4 2 

shown in Table 1. It is obvious if one graphs the results that as T increases 
the spectral density function tends toward a peak at 10 1/4 cycles per unit 
of time. At T = 4, the function is zero everywhere except at 10 1/4 cycles 
as it should be. Moreover, the peak value is 50 and this is the theoretical 
value. For T, any integral multiple of 4, the results will repeat. 

The values for the power spectral density at other frequencies when 
T < 4 are due to the model. Effectively the FFT calculates the Fourier 
series coefficients for a function which agrees with the given function from 
0 to T, and then repeats itself with period T . Clearly, in the given 
case, the two functions were not the same until T = 4. These values are 
usually referred to as the results of leakage. 

Since nothing is known about the values of the spectral density 
function between calculated values , discretion must be used in interpreting 
a single analyzed record. If a frequency component falls between calculated 
values of the spectral density function, the calculated values on each side 
will compensate for the missed frequency. The Fourier coefficients decay 

as — — ^ — , where f is the actual frequency and f is the calculated 
a c 
frequency. 

Studying two or more of these cases with T < 4 should give a 
better understanding than one case with a fixed T . 

The second example chosen was a square wave with amplitude 64 
and at 32 cycles per second. The sample rate was chosen to be 2048 samples 
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TABLE 1. SPECTRAL ANALYSIS OF x(t) = 10 cos 20. 5ir t 



N = 8; T 

= 1/8 

N = 16; T 

= 1/4 

N = 32; 

r *= 1/2 

N = 64; 

= 1 

N = 128; 

t = 2 

N - 256; 

t = 4 



G(f) 

PSD 

G(f) 

PSD 

G(f) 

PSD 

G(f) 

PSD 



G(f) 

PSD 

Frequency 

(Unit) 2 

Unit 2 /Hz 

(Unit) 2 

Unit 2 /Hz 

(Unit) 2 

Unit 2 /Hz 

(Unit) 2 

Unit 2 /Hz 

G(f) 

PSD 

(Unit) 2 

Unit 2 /Hz 

0.000 

0.865 

6.924 

0. 0737 

0.2948 

0.0607 

0.1215 

0.0968 

0.0968 

0. 0244 

0. 0122 


0.0000 

0.250 

0.500 

0.750 









0. 0246 

0.0123 




1.000 

1.250 







0. 0987 

0.0987 

0.0255 

0.0127 




1.500 

1.750 









0.0269 

0.0134 




2.000 





0.0655 

0.1310 

0.1046 

0. 1046 

0.0290 

0. 0145 




2.250 

2.500 

2.750 









0.0320 

0.0160 




3.000 







0.1158 

0.1158 

0. 0359 

0.0180 




3.250 

3.500 

3.750 









0.0411 

0. 0206 




4. 000 
4. 250 



0.2333 

0. 9332 

0.0837 

0.1674 

0.1346 

0. 1346 

0.0480 

0. 0240 




4.500 

4.750 









0.0571 

0.0286 




5.000 







0.1665 

0.1665 

0.0694 

0. 0347 




5.250 

5.500 

5.750 









0.0863 

0.0431 




6.000 

6.250 





0.1373 

0.2748 

0.2237 

0.2237 

0.1101 

0.0551 




6.500 

6.750 









0. 1448 

0. 0724 




7. 000 
7.250 







0. 3393 

0.3393 

0.1978 

0.0989 




7. 500 
7.750 









0. 2839 

0.1420 




8.000 

8.250 

5.881 

47. 049 

3.0507 

12.203 

0.3812 

0. 7625 

0.6322 

0.6322 

0. 4365 

0.2183 




8.500 

8.750 









0. 7431 

0.3715 




9. 000 







1.839 

1.839 

1.5010 

0. 7505 




9.250 
9. 500 
9. 750 









4.2989 

2.149 




10.000 





24. 4255 

48.8511 

41.541 

41.541 

39. 9042 

19.952 




10.250 











200.0 

50. 0000 

10.500 









41.1640 

20.582 


0.0000 

10.750 














11.000 







4.1890 

4. 1890 

4.7187 

2.359 




11.250 

11.500 

11.750 









1.7526 

0. 8763 




12. 00C 
12.250 



7. 4488 

29.7 

0. 3986 

0.7972 

0.7014 

0.7014 

0.9226 

0.4613 




12.500 









0. 5758 

0.2879 




12. 750 












0.0000 
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TABLE i. (Concluded) 



N - 8; T - 1/8 

N = 16 

T = 1/4 

N = 32; T = 1/2 

N = 64; t = 1 

N = 128; t = 2 

N = 256; t = 

1 


G(f) 

PSD 

G(f) 

PSD 

G(f) 

PSD 

G(f) 

PSD 



G(f) 

PSD 

Frequency 

(Unit) 2 

Unit 2 /Hz 

(Unit) 2 

Unit 2 /Hz 

(Unit) 2 

Unit 2 /Hz 

( Unit) 2 

Unit 2 /Hz 

G(f) 

PSD 

(Unit) 2 

Unit 2 /Hz 

13.000 







0,2600 

0.2600 

0.3977 

0.1988 


0.0000 

13.250 

13.500 

13.750 









0.2938 

0.1469 




14. 000 
14.250 





0.0698 

0.1396 

0.1285 

0.1285 

0.2276 

0.1138 




14.500 
14. 750 



- 






0.1828 

0.0914 




15.000 

15.250 







0.0739 

0.0739 

0.1510 

0.0755 




15.500 

15.750 









0. 1275 

0.0638 




16.000 

16.250 

0. 4995 

3.996 

0. 9435 

3.774 

0.0238 

0. 0477 

0. 0467 

0.0467 

0.1097 

0.0548 




16.500 

16.750 









0.0958 

0.0479 




17.000 







0. 0315 

0.0315 

0.0847 

0. 0424 




17.250 
17.500 
17. 750 









0.0758 

0.0379 




18.000 





0.0105 

0.0210 

0. 0223 

0.0223 

0.0684 

0.0342 




18.250 

18.500 

18.750 









0.0623 

0.0312 




19.000 







0.0164 

0.0164 

0.0572 

0.0286 




19.250 
19. 500 
19.750 









0.0528 

0.0264 




20.000 

20.250 



0.4362 

1.7446 

0.0052 

0.0105 

0.0124 

0.0124 

0. 0491 

0.0245 




20.500 

20.750 









0.0459 

0.0229 




21.000 

21.250 


- 





0. 0096 

0.0096 

0.0431 

0.0216 




21.500 

21.750 









0.0407 

0.0204 




22.000 





0.0028 

0. 0056 

0.0077 

0.0077 

0.0386 

0.0193 




22.250 
22.500 
22. 750 









0.0367 

0.0184 




23. 000 
23.250 







0.0063 

0.0063 

j 

0.0350 

0.0175 




23.500 
23. 750 








i 

0. 0336 

0.0168 


: 


24. 000 

0.0724 

0.5794 

! 0.2915 

1.6604 

0.0015 

0. 0030 

! 0. 0052 

0.0052 


0.0162 




25.000 







0.0044 

0. 0044 


0.0151 




26.000 







; 0.0038 

0.0038 


0.0142 




27.000 







0.0033 

0.0033 


0.0135 




28.000 



0.237 

0. 9487 

; 0.0005 

0.0009 

0.0030 

0.0030 


0.0130 


0 . 

0000 

TOTALS 




50.85 


51.33 








15 




per second with T = 1 second. The objective was to compare the accuracy 
of the FFT with the theoretical answers and also to compare it with the 
results from RAVAN and STAN for a similar time series. 

The calculated values of the PSD and the theoretical values for the 
first 10 components are tabulated in Table 2. The percentage error based 
on the maximum theoretical value is roughly 0. 04 percent. The actual 
value of the error decreases with higher frequency. 

A direct comparison with RAVAN and STAN values was not made 
since the example recorded in STAN 2 used a square wave of 100 cycles per 
second and 8000 samples per second. However, the actual errors in these 
cases increased with increasing frequency. With the FFT, the percentage 
error fell between the values of RAVAN and STAN for the lower frequencies 
For higher frequencies, the FFT answers had smaller percentage errors 
than either RAVAN or STAN, even though their sample rate was higher per 
cycle. 


TABLE 2. SQUARE WAVE (Amplitude 64) 


Freq. 

FFT (PSD) 

Freq. 

Theory (PSD) 

32 

3318. 76 

32 

3320. 0925 

96 

367.577 

96 

368. 899 

160 

131.577 

160 

132.8037 

224 

66. 487 

224 

67. 757 

288 

39. 763 

288 

40. 989 

352 

26.268 

352 

27.439 

416 

18.544 

416 

19.646 

480 

13.739 

480 

14, 756 

544 

10.572 

544 

11.488 

608 

8.400 

608 

9. 197 


2. ibid. 
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EXPERIMENTAL RESULTS 


Two channels of data from test AS-506 were used for several 
different studies. The first four of these studies were to show the effect 
of the sampling rate and the part of the records used. All four were taken 
from the two channels of data with a time slice from 60 seconds to 62. 048 
seconds. 

The first analysis used the first 1. 024 seconds of data from each 
channel with a sampling rate of 1000 samples per second. The predominant 
(largest PSD) frequencies and their magnitudes are shown in Table 3 for 
both channels. 

The second analysis used a time slice of i second from 60. 0 to 61. 0 
seconds of each channel. The sample rate was 1024 samples per second. 
The results of this analysis are shown in Table 4. 

Even though the sample rate and time interval were very nearly the 
same, some changes in the predominant frequencies did occur. No rigorous 
explanation of these differences is offered. Since the different lengths of 
record forced the calculation of the FFT at different discrete values of 
frequency, these variations may result from the sharpness of the peaks on 
the curves. 

The third and fourth analyses were made to further check the effect 
of sampling rate and length of record and also to check whether these 
particular time series were self-stationary with the given lengths of record. 

If the time series were self-stationary with time T , the results 
from any two segments of the record with length T would be the same. 

Table 5 shows the results for a time slice from 61. 024 to 62. 048 
seconds with a sample rate of 1000 per second. Table 6 shows the results 
of a time slice from 61. 0 to 62. 0 seconds with a sample rate of 1024 per 
second. 


Again, there is considerable variation between Tables 5 and 6 
because of the sample rates and lengths of record. Also, Tables 3 and 5 
differ as do Tables 4 and 6. This indicates that the time series are not 
self- stationary with records of length 1 second or 1. 024 seconds. It should 
be noted that the predominant frequency is basically the same for the given 
channel in all four cases. Many of the others occur in each case but their 
rankings are different. 
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TABLE 3, FFT — TIME FROM 60. 0 TO 61. 024 


Channel 1 

Channel 2 

Freq. 

Amp. (PSD) 

Freq. 

Amp. (PSD) 

186.523 

5.513 

187.500 

26.783 

156.250 

2.690 

281.250 

4.537 

177. 734 

2.172 

186.523 

4.263 

150.391 

2.094 

337.891 

3.973 

145.508 

2.042 

0. 977 

3.915 

151.367 

1.806 

213.867 

3.834 

168. 945 

1. 666 

235.352 

3.065 

187.500 

1.513 

182.617 

2.979 

162. 109 

1.474 

228.516 

2.675 

159. 180 

1.419 

207.031 

2.662 


TABLE 4. FFT - TIME FROM 60. 0 TO 61. 0 


Channel 1 

Channel 2 

Freq. 

Amp. (PSD) 

Freq. 

Amp. (PSD) 

187. 00 

8.015 

188.00 

28.641 

153.00 

2.729 

282.00 

5.52 

156.00 

2.262 

1.00 

5,157 

202.00 

1.831 

2.00 

4. 924 

152.00 

1.726 

344. 00 

4. 173 

178.00 

1.672 

226.00 

3.643 

163.00 

1.662 

215.00 

3.622 

6.00 

1.4914 

219.00 

3.315 

1.00 

1.424 

183.00 

2,939 

169.00 

1.135 

365.00 

1 

2.632 
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TABLE 5. FFT - TIME FROM 61. 024 TO 62. 048 


Channel 1 

Channel 2 

Freq, 

Amp. (PSD) 

Freq. 

Amp. (PSD) 

186.523 

10.208 

187.50 

21.869 

0.977 

6. 998 

188.477 

8.372 

158.203 

3.483 

2. 930 

6.717 

170. 898 

2.668 

215.820 

4.602 

163.086 

1.768 

226.562 

3.817 

10.953 

1.734 

204.102 

3.774 

147.461 

1.367 

220.703 

3.541 

180.664 

1.312 

239.258 

3.428 

195.312 

1.258 

214. 844 

3.378 

152.344 

1. 150 

347.656 

3.269 


TABLE 6. FFT - TIME FROM 61. 0 TO 62. 0 


Channel 1 

Channel 2 

Freq. 

Amp. (PSD) 

Freq. 

Amp. (PSD) 

187.00 

5. 838 

188.00 

23. 727 

186.00 

4.601 

189.00 

15.633 

160.00 

1.965 

215.00 

8.054 

194.00 

1.443 

282.00 

5.739 

7.00 

1.339 

281.00 

4.687 

159.00 

1.286 

233.00 

4; 063 

154. 00 

1.285 

249.00 

3.714 

60.00 

1.278 

205.00 

3.655 

180.00 

1.263 

172.00 

3.579 

140. 00 

1.235 

376. 00 

3.450 
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The next two analyses were made to compare the results of runs of 
length 2T with those of length T and also to compare the results of the 
same analyses performed by STAN. 

Table 7 shows the results of the analysis of 2. 048 seconds of data 
with a sampling rate of 1000 per second and Table 8 shows the results of 
analysing a 2 second record with sampling rate of 1024 samples per second. 
The agreement between these two cases seems to be better than for the 
shorter records. 

Tables 9 and 10 show the results of analysing the same data as 
Tables 7 and 8 but using STAN. The STAN program was set up to find the 
predominant frequencies up to 200 cycles per second. There is very good 
agreement between these methods as to the predominant frequency in each 
case. It should be noted that the slightly different sampling rate and length 
of records affected the results of the STAN program fully as much as it did 
the FFT. 

These results are also shown in Figures 2 through 9. Figure 2 
should be compared with Table 7, channel 1, and Figure 3 to Table 7, 
channel 2. Figure 4 corresponds to Table 8, channel 2. Figure 6 
corresponds to Table 9, channel 1, and Figure 7 to Table 9, channel 2. 
Figures 8 and 9 correspond to Table 10, channels 1 and 2, respectively. 

Two other experiments were performed to study the effect of the 
number of bits. When the data was changed from 10 bits to 6 bits, the 
change in the larger peaks was only about iO percent. Changing from 10 
bits to 2 bits gave about the same results for the major peak. These results 
are shown in Figures 10 through 13. From this, it appears that the number 
of bits used is not too significant. 

Two further analyses were performed. The purpose of these two 
was to check the results of the FFT when used to calculate the frequency 

response function of a system. Two time series, 100 e ^ and 
_2t 

100 e sin 12 7rt , were chosen. Since these functions are obviously not 
stationary, the spectral density function was modified. In this study, the 
quantity 2 X (f) X(f) was used and compared with the theoretical value of 
2 ! H (f ) 1 2 where H(f) is the frequency response corresponding to the given 
time response. The factor of 2 was included in order to use only positive 

frequencies. These results are given in Table 11 for 100 e and Table 
_2t 

12 for 100 e sin 127rt. The totals shown in the tables are for the integrals 
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TABLE 7. FFT — TIME FROM 60. 0 TO 62. 048 


Channel 1 

Channel 2 

Freq. 

Amp. (PSD) 

Freq. 

Amp.. (PSD) 

186.523 

6. 961 

187.500 

13.431 

0.488 

4.135 

187. 958 

7.586 

169.434 

1.474 

188.477 

4. 009 

1.465 

1.406 

188.956 

3.576 

0.977 

1.402 

215.332 

2. 904 

155. 762 

1.294 

0.488 

2.820 

152. 832 

1.181 

281.250 

2.675 

187.012 

1.1001 

223.145 

2. 549 

4.395 

1.063 

216.306 

2.064 

60.059 

0.967 

209.961 

1.997 


TABLE 8. FFT - TIME FROM 60. 0 TO 62. 0 


Channel 1 

Channel 2 

Freq. 

Amp. (PSD) 

Freq. 

Amp. (PSD) 

187.00 

3.764 

188.0 

13.898 

186.5 

3.457 

187.5 

7.615 

0.500 

2.023 

189.0 

6.251 

187.5 

1.635 

0.50 

5.301 

60.0 

1.568 

281.5 

4. 505 

156.0 

1.351 

188.5 

3. 162 

169.5 

1.140 

189.5 

2.352 

177.5 

1.049 

212.5 

2.236 

152.0 

1.034 

216.5 

2.230 

194.0 

0.944 

206.0 

2.156 
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TABLE 9. STAN - TIME FROM 60. 0 TO 62. 0 


Channel i 

Channel 2 

Freq. 

Amp. (PSD) 

Freq. 

Amp. (PSD) 

187.0 

10.005 

188.0 

30.393 

156. 0 

2.434 

182.0 

3.122 

170.0 

2.105 

196.0 

2.621 

6.0 

2.000 

193.0 

2.313 

153.0 

1.978 

165.0 

2.042 

163.0 

1.841 

8.0 

1.825 

159.0 

1.772 

172.0 

1.628 

178.0 

1.595 

176.0 

1.595 

60.0 

1.533 

60.0 

1.357 

195.0 

1.343 

155.0 

1.350 

1 

rABLE 10. STAN - ' 

riME FROM 60. 0 TO 62 . 048 

Channel 1 

Channel 2 

Freq. 

Amp. (PSD) 

Freq. 

Amp. (PSD) 

186.8 

8.898 

187.8 

26.646 

155.8 

1.753 

3.0 

3.240 

152.3 

1.555 

1.0 

2.615 

169.3 

1.407 

182.3 

2.235 

177.3 

1.399 

180.8 

1.820 

170.8 

1.369 

192.3 

1.713 

4.5 

1.313 

195.3 

1.711 

163.3 

1.151 

193.8 

1.402 

158.3 

1.036 

184.8 

1.266 

60.1 

1,022 

171.8 

1.203 
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TABLE ii. TOTAL POWER FOR 100 e 


Freq. 

1000 Samples/Sec 
1. 024 sec 

Theory 

0 

104. 17 

100.00 

0.977 

146.68 

145.30 

1. 953 

80. 57 

79. 81 

2.930 

46. 01 

45.58 

3. 906 

28.75 

28.47 

4. 883 

19.39 

19.21 

5.859 

13.87 

13.74 

6. 836 

10.38 

10.28 

7. 813 

8. 05 

7.97 

8.789 

6.41 

6.35 

9.766 

5.23 

5. 17 

Total 

497. 86 

500. 00 


of the area under the curves 2 X(f) X(f) and 2| H (f) 1 2 , respectively. 

It should be noted that the error in 2 |H(f) | 2 is 4 |H f | times as large 
as the error in | H (f) | . 

Hence, for the value of |H(f) | at 6 cycles per second in Table 12, 
the error should be approximately 1 part in 50 for 2 seconds of data and 
about 1 part in 500 000 when 16 seconds of data were used. 
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~2t 

TABLE 12. TOTAL POWER FOR 100 e sin 12vrt 


Freq. 

512 

Samples/ 
Sec 
2 sec 

512 

Samples/ 
Sec 
4 sec 

256 

Samples/ 
Sec 
4 sec 

512 

Samples/ 
Sec 
16 sec 

Theory 

0 

6.7392 

6.9865 

6. 9675 

6. 9865 

6. 9968 

1 

14.2477 

14.7744 

14.7353 

14. 7843 

14.7972 

2 

17.0152 

17. 6442 

17.6016 

17.6561 

17.6701 

3 

23.7952 

24. 6747 

24. 6245 

24.6914 

24.7082 

4 

42.7892 

44.3709 

44. 3037 

44.4006 

44.4642 

5 

131.7444 

136. 6147 

136.5003 

136.7063 

136.7443 

6 

1203.7800 

1248.2800 

1248.2704 

1249. 1187 

1249.1188 

7 

94.4052 

97. 8950 

97. 9904 

97. 9607 

97. 9291 

8 

21. 8684 

22.6768 

22. 7245 

22. 6920 

22.6760 

9 

8.5886 

8. 9061 

8. 9363 

8. 9121 

8.9020 

10 

4.2691 

4.4269 

4.4484 

4.4299 

4.4228 

11 

2.4271 

2.5169 

2. 5330 

2. 5186 

2.5132 

12 

1.5062 

1.5619 

1. 5747 

1.5630 

1.5587 

Total 






Power 

1242.70 

1244.74 

1244. 74 

i 

1246. 03 

1246.49 
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Figure 2. FFT — Power spectral density, Channel 1, 1000 samples per second 
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Figure 4. FFT — Power spectral density, Channel i, 1024 samples per second 
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Figuer 5. FFT — Power spectral density, Channel 2, 1024 samples per second 
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Figure 6. STAN — Power spectral density 



Channel 1, 1000 samples per second 
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Figure 7. STAN — Power spectral density, Channel 2, 1000 samples per second. 
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Figure 8. STAN — Power spectral density, Channel 1, i024 samples per second. 
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Figure 9. STAN — Power spectral density, Channel 2, 1024 samples per second. 
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Figure li. Power spectral density from six-bit data, Channel 2. 






The results for these two cases are very accurate. The maximum 
percentage error for |H(f)| is for the 0 frequency value of H(f) in 
Table 1 1 and this is less than 2 percent. The maximum percentage error 
relative to the maximum value of |H(f) | for the other terms in Table ii 
is less than 0. 5 percent. In Table 12, the maximum percentage error in 
I H (f) | relative to the maximum value of | H (f ) | is less than 0. 2 percent. 


INTERPRETATION 


The spectral density function is calculated at discrete points. In 
general, nothing is known about its value between these points. The graphs 
in Figures 2 through 9 are the products of spectral density function and 
bandwidth. 

Without prior knowledge of the spectral density function, one cannot 
tell whether indicated energies at adjacent points result from one or more 
frequency components. In this case, a more detailed analysis is required. 
Frequently, the analysis of a longer segment of the time series will suffice. 
This is verified by Table 1. 

In the cases studied here, it appears that most of the smooth peaks 
resulted from one frequency falling between the calculated values. In the 
two cases of actual observed data, it was obvious that small changes in the 
points where the FFT was calculated made large changes in the maximum 
values calculated. In other words, the maximum calculated value may not 
be a very good approximation of the true maximum value in the neighborhood 
of that point. However, with the recommended sampling rates and length 
of series analyzed, the values of the spectral density function are as good 
approximations of the true values as any other discrete method. 


TYPICAL HARDWARE AND SOFTWARE CHARACTERISTICS 

The time required to calculate the FFT and the power spectral 
density depends on the number of values of the time series, the number of 
bits used to express the value of the function, and the computer and display 
used. 
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The number of computation operations required to perform the FFT 
and N data points is proportional to N log 2 N versus N 2 for the 
classical Blackman-Tuckey approach. When N is greater than 1023, the 
time can be about 100 times faster for the FFT approach then for the 
classical approach. An even greater time savings can be achieved utilizing 
hardware to calculate the FFT, which is called the Fast Fourier Analyzer 
(FFA) . 


Special purpose hardware, such as the FFA hardware, can be 
attached to the mainframe computer as a peripheral to one of the mainframe 
input/output (I/O) channels. The mainframe computer controls the FFA 
hardware via priority interrupts or other program^ logic. The FFA hard- 
ware has a separate core memory to allow simultaneous operations to be 
performed by the mainframe computer and FFA hardware. This is a very 
important advantage for many real-time data processing applications such 
as image enchancement, spectral analysis, radar, sonar, and vibration 
analysis. See summary Table 15 for some typical FFA's. 

Since the FFA typically produces a squared quantity summed over 
a large data population N , the number of bits to resolve the correct 
answer is defined by the following equation: 


R = J + 2(K) + 1 


where J is the number of bits that will define the maximum population N 
to be considered by the FFA. K is the A/D bit word size not including the 
sign. Therefore, the FFA hardware specifications that define the complex 
arithmetic bit register resolution to maintain accuracy without data 
compression is R bits. 

Often, the test engineer wants to know how many data channels may 
be analyzed by the FFA in real time. Since this will depend on the specific 
FFA and computer used, these timing equations are given as a guide to aid 
the engineer in this evaluation. 


INPUT 

T. 

l 

= HBN + P. 

jit sec 


X 

OUTPUT 

T 

= HBN + P 

ju sec 


o 

s 

FFT 

T 

f 

= 4N log N + 2N + P 

s 

j u sec 
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MULT 


H sec 


Spectral 

Density 


T = 6N + P 
m r 



T +T 


m 


N is the number of data points (must be a power of two) of the 
real or imaginary data array; 

B is the number of bytes in each element of the data array; 

P. is the program CPU overhead to acquire input data, convert 

to engineering unit, and load FFA data memory; 

P is the CPU program overhead to service the priority interrupt 
service subroutine, store complex results, and decrement 
counters and other program logic to orderly proceed through 
the FFT; 

H is the channel transfer time per byte; 

P is the CPU program overhead to convert the FFT real and 

imaginary array to power spectral densities, and conversion 
for output display. 

Tables 13 and 14 show a typical test case using two different 
computers to perform the FFT with and without FFA. To perform on-line 
quick-look analysis, the most time consuming part is the integration time 
of the FFA (Table 14). In our investigation, the next time consuming item 
was the display time of the CRT. This limits the number of the data 
channels that may be analyzed during real time to one data channel every 
13 seconds for the SDS-930 computer system and one channel for every 
five seconds for the SIGMA- 5 computer; real time could not be considered 
for most vibration applications because the software speed of the FFT would 
be too slow. Without the FFA, the time required for analysis and display 
for one channel is 98 seconds with the SDS-930 computer and 12. 2 seconds 
for the SIGMA -5 computer. The type of mainframe computer used will 
influence the choice of the FFA hardware. Also, a fast and more advanced 
computer may not require FFA hardware for some applications. 
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TABLE 13. TYPICAL CHARACTERISTICS OF ENVIRONMENTAL TESTS 


Type 

Of 

Test 

Typical Time 
Duration/Test 
(sec) 

Highest 

Frequency 

Content/Ch. 

(Hz) 

Typical Sampling 
Rate Channel 

No Channels 
Required 
Total/On-Line 

System Recording/ 
Thoroughput Rate 
( On-Line) a 

Sine 

0 - 540 

2000 



14/2 

115 200 

Sweep 

540 

1000 



28/12 

115 200 


540 

500 


■ 

32/12 

80 000 


540 

30 



32/12 

4 800 

Sine 

0-10 

3000 


10 000 

14/2 

115 200 

Dwell 

10 

1000 


5 000 

28/12 

115 200 


10 

500 


2 500 

32/12 

80 000 


10 

30 


150 

32/12 

4 800 

Random 

0-10 

5000 

10 000 

12 000 

11/2 

115 200 


10 

1 

2000 

4 000 

6 000 

28/2 

115 200 


a. 


Data samples/second 














TABLE 14. FFT TIMING ESTIMATES 


Operation 

SDS-930 

XDS-SIGMA 5 

N = 2048 
(sec) 

N = 4096 
(sec) 

N = 2048 
(sec) 

N = 4096 
(sec) 

1 

Acquire x(t) in CPU 
mainframe memory 

■ 

0. 5 

1.0 

0. 1 

0.2 

Transform to 
engineering units 

T. 

l 

0.25 

0. 5 

0. 1 

0. 2 

Transform to 
zero mean 

■ 

0. 1 

0. 2 

0. 05 

0. 1 

Spectra Computation 
by software without 

FFA a 

T 

f 

s 

90 

270 

10 

31 1 

b 

Load FFA memory 

fl 

0.006 

0. 012 

0.003 

0. 006 

FFA solution time (FFT 
service subroutine plus 
overhead and spectra 

conversion) ^ 

P 

s 

5 

15 

2. 5 

6.5 

CRT display overhead 

P 

r 

7 

14 

2 

4 

Total time with FFA 
simultaneous input/ 
output throughput^ 


12 

27. 77 

4.753 

11. 0006 

Total time 
simultaneous input/ 
output with software 

FFT a 


97. 856 

285. 712 

12. 253 

35. 506 


a. Software option 

b. Hardware option 
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PROBLEM AREAS 


Some work was done to determine the effect of the number of bits 
used in the quantization of the data on the calculated values of the FFT and 
the power spectral density function. The few cases studied seem to indicate 
that the number of bits used was not too significant. Some improvement in 
the speed of calculation usually follows when fewer bits are used. However, 
scaling problems are introduced and the effect of the small changes in FFT 
and power spectra may cause large changes in the autocorrelation function 
and other related functions. This area needs more study. 

Another area of concern is the autocorrelation and cross-correlation 
functions. Most of the literature recommends some type of filter of the 
original data. With deterministic data, the results were good without 
these filters. Study of the effect of these filters on other types of data is 
indicated. 

The transfer function and phase angle need to be investigated further. 
The use of these functions without separation of frequencies seems to give 
almost meaningless results. 

Figures 14 and 15 show the raw data that were analyzed. 


CONCLUSION 


The FFT that has been defined and used in this report does yield 
good estimates in the frequency domain. These estimates are given at 
discrete intervals (Af). The FFT, as all other digital methods such as 
STAN or RAVAN, does not define energy between two adjacent frequency 
intervals. Values between adjacent points can be interpolated; however, 
caution has to be taken in the interpretation of these interpolated results. 
When more resolution is desired, a longer data record can generally be 
used to obtain this resolution. The FFT does yield results at least equal 
to previously used techniques such as RAVAN and STAN. More important, 
the speed advantage of the FFT is several magnitudes greater than the 
Blackman-Tuckey techniques used by RAVAN. The most important feature 
of the FFT is that computational algorithm hardware and core memory 
(FFA) can be added to existing computers to perform rapid real-time 
signal or data processing. In our investigation, the number of bits that 
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Figure 15. Raw data, Channel 2 



defined each data sample was not as critical as the sampling rate. Figures 10 
and 11 showed that reducing the A/D size from 11 bits to 6 bits only reduced 
the power spectral density (PSD) magnitude of the major five predominant 
frequencies by 10 percent. When the A/D size was reduced from 11 bits to 
2 bits, Figures 12 and 13, the PSD magnitude of the predominant frequency 
was degraded about 10 percent. As expected, the number of predominant 
frequencies that are detectable are fewer and the lower amplitudes are 
degraded the most when two bits are used. 

From the foregoing studies it can be concluded that quick-look 
vibrational analysis of a few channels with a digital computer is feasible for 
structural testing. A hardware Fast Fourier Analyzer (FFA) connected to 
the computer as a peripheral reduces considerably the analysis time, 
particularly in case a relatively slow computer is used. It is expected that 
with a modern medium size digital computer with an FFA, a complete 
quick-look vibrational analysis of a signal can be computed and displayed 
every 5 seconds. If one assumes that a waiting time of 1 minute is 
acceptable, 12 channels could be monitored simultaneously for vibrational 
analysis of structures. The study also showed that the time to display the 
data on the CRT must not be neglected, and that proper selection of the 
display unit has to be considered. 

The FFT was used to calculate the total power in two response 

functions (time function). These two functions were 100 e 10t and 
-2t 

100 e sin 12irt. These were chosen to represent two realizable systems. 
The calculated values and the theoretical values were in close agreement. 

The results are summarized in Tables 11 and 12. 
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APPENDIX 


FAST FOURIER TRANSFORM 


The spectral density functions can be defined directly by taking Fourier 
transforms [1]. There are certain theoretical requirements necessary for the 
existence of these transforms; they will be assumed to be met. For those 
interested,, there are many rigorous treatments of the Fourier transform. The 
spectral decomposition of a time series x(t) will be developed by assuming that 
it has a complex Fourier transform X(f) such that 


X(f) = / x(t) e~i 2n ft dt 

-OO 

and conversely 


x(t) 


/ X(f) e j27rft df 


A sufficient condition for the existence of these integrals is that x(t) and its 
derivative x(t) be piecewise continuous in every finite interval (a,b) and the 
x(t) be integrable on (-<», °°) . These conditions can be satisfied (and usually 
are) in practical problems by setting x(t) equal to zero outside some fixed range 
of t [that is, making x(t) a finite time series] . 

Similarly, if there exists a second time function y(t) , it will be assumed 
to have a complex Fourier transform Y(f) such that 

Y(f) = / y(t) e -j2ir ft dt 

— OO 

and 

y(t) = / Y(f) e j27rft df 
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In practice, x(t) and y(t) are assumed to be zero outside some range of t 
and that there are no frequencies above some finite frequency F, Hence, the 
integrals become 

X(f) = f x(t) e~ j2w ft df 
0 

x(t) = / X(f)e j27r ft dt 

-F 

Y(f) = f y(t) e~j 2w ft dt 
0 

and 

F 

y(t) = f Y(f) e i27rft df 

-F 

Also, in practical problems the time series x(t)' must be sampled at a certain 
increment and the sampled values used. This requires the integrals to be 
redefined as sums of products. Usually, equal spacing is used for time and 
X(f) and Y(f) are calculated at equally spaced points of frequency. Using 

the trigonometric form of e J , these integrals become 



N-i 



X(f.) = (At) 

2 

K=0 

X( V 

[cos 27 r f.t. - j sin 2tt f.t, ) 
i k J i k' 


N-l 



x(y = (Af) 

2 

i=0 

X(f.) 

[ cos 27T Tt^ + j sin 2tt f.t^) 

and similarly for Y ( f) and 

y(t) . 



If X(f) and Y(f) are calculated for a discrete set of t values, they 
are referred to as Discrete Fourier Transforms or DFT. When the time series 
x(t) is a discrete set of equally spaced numbers, there is a process of 
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T 

multiplication that permits faster calculation of X(f) . If tj = — where 


N = 2 for some integer and 

-j27 r fiti . , , 

e is denoted by a , 

X(f) can be written as 


4f = ^, 2n f,t, = 2* (1) (|) = f . Then if 

f^i ^i 

then e = a . Using this notation 


X(0) 


1 1 


x(t) 



2 N-i 



X(Af) 


i a a ... a 


x(t.) 


= 

< 2 4 2(N-i) 

1 a a * , . a 


9 



' J 2 



X[ (N-l) Af] 


1 a (M ' 1) . a (N - 1) 


j; (t N-iL 


T T 

or X (f) = Ax (t), where each value of X(f) must be multiplied by (At) 
for the true value. 


If this method is used with the following arithmetic, it is called the Fast 
Fourier Transform or FFT. The speed of calculation is obtained by factoring 
A in such a way as to minimize the complex multiplication and addition. It 

_27T 0 

should be remembered that e = e = i for any integer k and since 


_ i^i 

a = e ^ can be represented as a * where t is the remainder when 

kl is divided by N. Using this notation, factoring for N = 4 and N = 8 is 
given. If N = 4 (i.e., x 0 , Xj, x 2 , x 3 form the discrete time series) , the 
matrix notation for the Discrete Fourier Transform is 


where 


X(f) = At 



"l l 1 i ' 


x 0 


X(0) 

1 a a 2 a 3 


*l 


X(Af) 

1 a 2 a 4 a 6 


x 2 


X ( 2 A f ) 

1 a 3 a 6 a 9 


x 3 

! 

X(3Af) 


i7r 4k + 1 1 . , . . „ 

= — — , a = a and 1 = 0, i, 2, 3. 
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The Fast Fourier Transform factors A as follows: 


1 

a 4 

0 

0 



1 

0 

1 

0 

1 

a 2 

0 

0 



0 

i 

0 

1 

0 

0 

1 

a 4 

and 

R = 

1 

0 

Ol 2 

0 

0 

0 

1 

a 2 



0 

a 

0 

a 3 


but 


1 

1 

i 

1 

1 

a 2 

Ol 4 

a 6 

1 

Ol 

a 2 

a 3 

i 

a 3 

Ol 3 

a 9 


The two middle rows have to be interchanged. Hence, if one does the matrix 
T T T 

multiplication y (t) =R[x(t)], then y (t) = Qy (t) . The rows must be 
decoded. Using binary subscripts and starting with 0, 


Roo 

"-*■ Roo 

Roi 

Rio 

R i0 

Roi 

Rn 

— Rjj 


In other words the second element of Z would be the third element of AX 

T 

and the third element of Z would be the second element of AX . The first 
and fourth (zero and third) would be in the same locations. 

If N = 8 (i.e. , x 0 , x ls x 2 , x 3 , x 4 , x 5 , x e , x form the discrete time 
series) , the matrix multiplication for the Discrete Fourier Transform is 
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and X[ (k-i)Af] is the vector dot product of the kth row of the matrix and the 
column of x values. The powers of a are not reduced module 8 for reasons 

( 

that will became apparent \a = e ) . This same matrix is used for the 
Fast Fourier Transform but is factored into three matrices and multiplied in 
that way. 

The matrix factorization is done in the following way. Take the first 
two rows (2x4 = 8) and divide each one into two equal parts. The first 
half of the first row is the diagonal of a 4x4 submatrix in the upper left 
corner. The second half of the first row is the diagonal of a 4x4 matrix 
in the upper right corner. Then the first half of the second row is the diagonal 
of a 4 x 4 matrix in the lower left corner and the remainder of the second row is 
the diagonal of a 4x 4 submatrix in the lower right. In short, one of the factors 
is 

1 0 0 0 i 0 00 

010 00 1 00 

001 00 0 10 

000 10 0 01 

100 0 a 4 0 00 

0 a 0 o 0 o; 5 0 0 

0 0 a 2 0 0 0 Q! 6 0 

000 a 3 0 0 0 a 7 
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Now, the other 

factor must be found. It is 
— • 





1 

1 

1 

1 

0 

0 

0 

0 


1 

a 1 

a 4 

a 6 

0 

0 

0 

0 


1 

a 4 

a 8 

a 12 

0 

0 

0 

0 


1 


a 12 

a 18 

0 

0 

1 

0 

Q = 

0 

0 

0 

0 

1 

i 

i 

1 


0 

0 

0 

0 

i 

a 2 

a 4 

a 6 


0 

0 

0 

0 

1 

a 4 

a 8 

a 12 


0 

0 

0 

0 

1 

a 6 

a 12 

a 18 


However QR interchanges rows of A so that they occur R 1( R 3 , R g , R 7 , R 2 , 
R 4 , R e , R 8 . (Perhaps one would rather say R 0 , R 2 , R 4 , Rg, R 1# R3, R g , R 7 , 
since if binary subscripts are used this is more consistent. ) 

Now, Q must be factored. If these are designated as Q1R1, 



i 

0 

1 

0 


0 

1 

0 

1 

0 

i 

0 

a 4 

0 

0 

a 2 

0 

a 6 


0 1 

B 


where B is the matrix in the upper left hand corner. Basically the same 
procedure as used to get R from A is used on each 4x4 matrix on the diag- 
onal of Q. Then, Q is a matrix with 2x2 matrices on the diagonal so that 


Qi - 


1 

1 

0 

0 

0 

0 

0 

0 


a 8 0 

a 4 0 

0 1 

0 i 

0 0 

0 0 

0 0 

0 0 


0 0 

0 0 

a 8 0 

a 4 0 

0 1 

0 1 

0 0 

0 0 


0 0 

0 0 

0 0 

0 0 

a 8 0 

a 4 0 

0 1 

0 1 


0 

0 


0 

0 

0 

0 


a 


a 


4 
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Actually, one should remember that a; 8 = i, a; 4 = -i and in general 

^ (8k + n) _ a n w h ere n = o, 1, . . . , 7; also (Q, R, ) interchanges rows of 
Q. The old rows of Q are now ordered R^, R g ', R 2 ', R 4 ', R 5 ’, R 7 ', R e ', R 8 ', 
where R ’ was the kth row in Q. Hence, the rows are now ordered quite dif- 

K. 

ferently from the original matrix A. In matrix notation, the relationship can 
be written A = PQR = PCPiQjRj) R, where P is a matrix obtained from the 
identity matrix by doing to it what must be done to R to give A; that is, first 
row is left alone, second row becomes third, third row becomes fifth, fourth 
row becomes seventh, fifth row becomes second, sixth row becomes fourth, 
and seventh row becomes sixth. This last row is unchanged. This gives 


P 


1 

0 

0 

0 

0 

0 

0 

0 


Now, similarly 


Pi 


1 

0 

0 

0 

0 

0 

0 

0 


0 

0 

i 

0 

0 

0 

0 

0 


0 

0 

1 

0 

0 

0 

0 

0 


0 0 
0 0 
0 0 
0 0 
1 0 
0 0 
0 1 
0 0 

0 0 
i o 
0 0 
0 i 
0 0 
0 0 
0 0 
0 0 


0 0 
i o 
0 0 
0 1 
0 0 
0 0 
0 0 
0 0 

0 0 
0 0 
0 0 
0 0 
1 0 
0 0 
0 i 
0 0 


0 0 
0 0 
0 0 
0 0 
0 0 
i 0 
0 0 

0 i 

0 0 
0 0 
0 0 
0 0 
0 0 

1 0 
0 0 
0 i 


and 
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PPJ = 


1 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

i 

0 

0 

0 


0 

0 

1 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

1 

0 


0 

1 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

i 

0 

0 


0 

0 

0 

i 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

1 


In other terms the first and eighth rows of R are the first and eighth 

rows of A, the fifth row of the product is the second' row of A, etc. An easier 
way to remember is to use binary subscripts on the rows and start with 0. Then 
they are R 000 , R„0i> R oio> R ou> R ioo> R ioi> R no> R iii- Reverse the order of the 
ones and zeros. The rows in the product Q^R^ R must be reordered in that 
manner; that is, 


In A 
^ R ooo 
R ooi 
R oio 
R oii 
R 100 
R 101 
R uo 

R lil 

For clarity the fourth row starting wi 
of A. 


In Product 

R ooo 

R 100 

R 0i0 

R liO 

- R 00i 

- R 10i 

- R on 

R iii 

the 0 row of QjRjR is the first row 


The order of multiplication then becomes 


R X 
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and 


T T 

U = Q t W 


T T 

X (f) = PPiU 

T T 

or X (f) is U decoded by the process explained above. The factor (At) 
is usually omitted until the last step. 

There are several advantages using this method. One of these is as 
follows. X(f) is permitted to be complex. If one has two real time series, 
let Z(t.) = x(t ) + j y(t.) . Find the FFT of z(t) or Z(f). Then 

Z(kAf) + Z[(N-k)Af] 

2 

Z(kAf) - Z[(N-k) Af] 

2j 

where Z is the complex conjugate of Z . 

Another special application is the following. If x(t) is a discrete 
(equally spaced) time series of length 2 , let 


X ± (t) = (X , x 2 , x 4 , . . , 

• X N-2> 

x 2 (t> = < X 1- x 3' V • • 

■ X N-1> 


be the two discrete time series of length 2 formed as shown. Then enter 
Xi(t) as the real part of x(t) and X 2 (t) as the imaginary part of x(t) in the 
matrix. Now, noting that the time spacing in xj(t) and x 2 (t) is 2At, the 
FFT of x(t) can be obtained by one pass through FFT process. Let A t (f) 
be the transform of x^t) and A 2 (f) be the transform of x 2 (t) . Then, for 
the original series x(t) , 


X(kAf) 

and 

Y(kAf) 
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X(kAf) 


j. 

2 


Aj(kAf) +A 2 (kAf)e 


ftn k 
N 


for k < — and 

X [ ( N-k) A f ] = X(kAf) 
for (N-k) > ^ . 

Since the correlation function is the inverse Fourier transform of the 
power spectral density function, the FFT routine may be used to calculate the 
correlation function. Actually, the entries in the matrix notation would have 
positive powers. However, the problem is easily solved by complex arithmetic. 

Most programs use e^ 2lT ft rather than e~^ 2ir ft as used in this report. If this 
is the case, then the actual Fourier transform is the complex conjugate of the 
calculated FFT. If this is the case, then using S (f) as the entry for x(t) 

X 

yields the autocorrelation function. In matrix notation, using C (t) for the 

XX 

correlation function, 


C (0) 

XX 


1 1 . . . 

i n 


X(0) 

C xx (t '> 


j 27T 

. + v 

1 e 

j47T j(N-l)27T 

n N 

e e 


X(Af) 

- 

= (Af) 

• 



* 



(N-l)j2 7T 

j(N--i) 2 2 tt 



C xx (t N-l ) 


< N 
1 e 

N 

e 


X[ (N-l)A f ] 



— 




or 

T T 

C (t) = BX (f) 

XX 

where each element of B is the complex conjugate of the matrix A in the 
calculation of the FFT. In most subroutines, the matrix B , given above is 
the one used. Logic steps easily change the calculated values into the correct 
Fourier transform. 
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